Datos satelitales en GRASS GIS

Los datos satelitales en general vienen en formato raster, por lo tanto aplican las mismas reglas.

Los comandos i.* se orientan explícitamente al procesamiento de datos satelitales aunque algunos puedan usarse para otros datos raster.

Nota

Para más detalles, ver el manual de Imagery Intro y la wiki sobre Image Processing.

Datos para esta sesión

Escenas Landsat 8 (OLI)

• Fechas: 14/01/2020 y 02/03/2020
• Path/Row: 229/082
• CRS: UTM zona 20 N (EPSG:32620)

Descargar las escenas L8 14/01/2020 (979Mb) y L8 02/03/2020 (880Mb) y moverlas a $HOME/gisdata/landsat_data. No descomprimir!

Historia de la mision Landsat

Lanzamientos de satélites Lansat desde 1972

El sistema de escáner multiespectral (MSS) a bordo del Landsats 1-5 disponía de cuatro bandas. El Thematic Mapper (TM) a bordo de Landsats 4 y 5 tenía siete bandas. El Enhanced Thematic Mapper Plus (ETM+) del Landsat 7 tiene 8 bandas y los Landsats 8 y 9 tienen 11 bandas. Fuente: https://landsat.gsfc.nasa.gov/satellites/landsat-9/landsat-9-bands/.

Comparación entre las bandas de todos los satélites Landsat

Nota

Más detalles sobre las misiones Landsat pueden encontrarse en: https://www.usgs.gov/landsat-missions

Manos a la obra

Iniciamos GRASS GIS

Iniciamos GRASS GIS en posgar2007_4_cba/PERMANENT

import os

# data directory
homedir = os.path.expanduser('~')

# GRASS GIS database variables
grassdata = os.path.join(homedir, "grassdata")
project = "posgar2007_4_cba"
mapset = "PERMANENT"

# import standard Python packages we need
import sys
import subprocess

# ask GRASS GIS where its Python packages are to be able to run it from the notebook
sys.path.append(
    subprocess.check_output(["grass", "--config", "python_path"], text=True).strip()
)

Ahora sí, estamos listos para importar los paquetes de GRASS e iniciar una sesión:

# import the GRASS GIS packages we need
import grass.script as gs
import grass.jupyter as gj

# Start the GRASS GIS Session
session = gj.init(grassdata, project, mapset)

Corroboramos la proyección

# check the CRS
gs.read_command("g.proj", flags="p")

Crear un nuevo mapset

Creamos un nuevo mapset llamado landsat8:

# Create a new mapset
gs.create_environment(grassdata, project, mapset="landsat8")

Listamos los mapsets accesibles

# list all the mapsets in the search path
gs.mapsets()

Listamos los mapas vectoriales disponibles

# list vector maps in all mapsets in the search path
gs.list_grouped(type="vector")

Región de interés

Extraemos el radio urbano de Córdoba

# extract Cordoba urban area from `radios_urbanos`
gs.run_command("v.extract", input="radios_urbanos", where="nombre == 'CORDOBA'", output="radio_urbano_cba")

Establecemos la región computacional al radio urbano de Córdoba

# set the computational region to the extent of Cordoba urban area
gs.run_command("g.region", flags="p", vector="radio_urbano_cba")

Descargar e importar los datos L8

Instalar la extensión i.landsat:

# install i.landsat toolset
g.extension extension=i.landsat

Buscar escenas de Landsat 8 disponibles

# search for Landsat 8 scenes
i.landsat.download -l settings=$HOME/gisdata/USGS_SETTING.txt \
  dataset=landsat_8_c1 clouds=35 \
  start='2019-10-27' end='2020-03-15'

NO EJECUTAR! Descargar las escenas seleccionadas

# download selected scenes
# i.landsat.download settings=$HOME/gisdata/USGS_SETTING.txt \
#   id=LC82290822020062LGN00,LC82290822020014LGN00 \
#   output=$HOME/gisdata/landsat_data

Imprimir las bandas dentro de la carpeta

# print all landsat bands within landsat_data folder
gs.run_command(i.landsat.import -p input=$HOME/gisdata/landsat_data)

Imprimir sólo las bandas seleccionadas con un patrón

# print a selection of bands - might be sloooow
gs.run_command(i.landsat.import -p input=$HOME/gisdata/landsat_data pattern='B(2|3|4|5|6|8)')

Importar bandas, recortar y reproyectar al vuelo

# import all bands, subset to region and reproject
gs.run_command(i.landsat.import -r input=$HOME/gisdata/landsat_data extent=region)

Listar bandas importadas y revisar metadatos

# list raster maps
gs.list_grouped(type=raster)["landsat8"]

# check metadata of some imported bands
gs.raster_info(map=LC08_L1TP_229082_20200114_20200127_01_T1_B4)
gs.raster_info(map=LC08_L1TP_229082_20200114_20200127_01_T1_B8)

Pre-procesamiento de datos satelitales

Workflow de pre-procesamiento de datos satelitales

De número digital (ND) a reflectancia y temperatura

• Los datos L8 OLI vienen en 16-bit con rango de datos entre 0 y 65535.
• i.landsat.toar convierte ND en reflectancia TOA (y temperatura de brillo) para todos los sensores Landsat. Opcionalmente proporciona reflectancia de superficie (BOA) después de la corrección DOS.
• i.atcorr proporciona un método de corrección atmosférica más complejo para gran variedad de sensores (S6).

Definir region computacional a banda de 30m

# set the region to a 30m band
gs.run_command("g.region", raster="LC08_L1TP_229082_20200114_20200127_01_T1_B4", flags="p")

Convertir DN a reflectancia superficial y temperatura - método DOS

# convert from DN to surface reflectance and temperature - requires to uncompress data locally
gs.run_command(i.landsat.toar input="LC08_L1TP_229082_20200114_20200127_01_T1_B"   output="LC08_229082_20200114_toar_B" sensor="oli8" metfile="$HOME/gisdata/landsat_data/LC08_L1TP_229082_20200114_20200127_01_T1_MTL.txt"
method="dos1")

Corroborar info antes y después de la conversión para una banda

# list output maps
g.list type=raster mapset=. pattern="*toar*"

# check info before and after for one band
r.info map=LC08_L1TP_229082_20200114_20200127_01_T1_B3
r.info map=LC08_229082_20200114_toar_B3

Banda 10 de L8 con la paleta de colores kelvin

Ahora, sigan los mismos pasos para la escena del 02/03/2020. ¿Qué notan de diferente?

Ajuste de color y composiciones RGB

Ajuste de colores para una composición RGB color natural

# enhance the colors
gs.run_command(i.colors.enhance red=LC08_229082_20200114_toar_B4   green=LC08_229082_20200114_toar_B3 blue=LC08_229082_20200114_toar_B2 strength=95)

Mostrar la combinación RGB

# display RGB
rgb = gj.

Seguir los mismos pasos para una composición falso color 543. Sobre qué bandas debieran realizar el ajuste?

Composición color natural 432

Composición falso color 543

Enmascarado de nubes con banda QA

• Landsat 8 proporciona una banda de calidad (quality assessment, QA) con valores enteros de 16 bits que representan las combinaciones de superficie, atmósfera y condiciones del sensor que pueden afectar la utilidad general de un determinado pixel.
• La extensión i.landsat.qa reclasifica la banda QA de Landsat 8 de acuerdo a la calidad del pixel.

Nota

Más información sobre la banda QA de L8 en la guía de usuario.

Crear las reglas para identificar las nubes y sombras de nubes

# create a rule set
gs.run_command(i.landsat.qa collection=1 cloud_shadow_confidence="Medium,High"   cloud_confidence="Medium,High" output=Cloud_Mask_rules.txt)

Reclasificar la banda QA en función de las reglas

# reclass the BQA band based on the rule set created
gs.run_command(r.reclass input=LC08_L1TP_229082_20200114_20200127_01_T1_BQA
  output=LC08_229082_20200114_Cloud_Mask rules=Cloud_Mask_rules.txt)

Reporte del porcentaje de nubes y sombras

# report % of clouds and shadows
print(gs.read_command(r.report -e map=LC08_229082_20200114_Cloud_Mask units=p))

Mostrar el mapa reclasificado

# display reclassified map over RGB
d.mon wx0
d.rgb \
  red=LC08_229082_20200114_toar_B4 \
  green=LC08_229082_20200114_toar_B3 \
  blue=LC08_229082_20200114_toar_B2
d.rast LC08_229082_20200114_Cloud_Mask

Comparar visualmente la cobertura de nubes con la composición RGB 543.

Composición falso color

Máscara de nubes

Fusión de datos/Pansharpening

Vamos a usar la banda PAN (15 m) para mejorar la definición de las bandas espectrales de 30 m, por medio de: i.fusion.hpf, que aplica un método de adición basado en un filtro de paso alto. Otros métodos están implementados en i.pansharpen.

Instalar la extensión i.fusion.hpf

# Install the reqquired addon
gs.run_command(g.extension extension=i.fusion.hpf)

Cambiar la región a la banda PAN

# Set the region to PAN band (15m)
gs.run_command(g.region -p raster=LC08_229082_20200114_toar_B8)

Ejecutar la fusión

# Apply the fusion based on high pass filter
gs.run_command(i.fusion.hpf -l -c pan=LC08_229082_20200114_toar_B8, msx=`g.list type=raster mapset=. pattern=*_toar_B[1-7] separator=,`, suffix=_hpf, center=high, modulation=max, trim=0.0)

Listar los mapas resultantes usando un patrón de búsqueda

# list the fused maps
gs.list_grouped(type=raster, pattern=*_hpf)["landsat8"]

Visualizar las diferencias

Datos originales 30 m y datos fusionados 15 m

# display original and fused maps

Índices de agua y vegetación

Establecer la máscara de nubes para evitar el cómputo sobre las nubes

# Set the cloud mask to avoid computing over clouds
gs.run_command(r.mask raster=LC08_229082_20200114_Cloud_Mask)

Calcular el NDVI y establecer la paleta de colores

# Compute NDVI
r.mapcalc \
  expression="LC08_229082_20200114_NDVI = \
  (LC08_229082_20200114_toar_B5_hpf - LC08_229082_20200114_toar_B4_hpf) / \
  (LC08_229082_20200114_toar_B5_hpf + LC08_229082_20200114_toar_B4_hpf) * 1.0"
# Set the color palette
r.colors map=LC08_229082_20200114_NDVI color=ndvi

Calcular NDWI y establecer la paleta de colores

# Compute NDWI
r.mapcalc expression="LC08_229082_20200114_NDWI = \
  (LC08_229082_20200114_toar_B5_hpf - LC08_229082_20200114_toar_B6_hpf) / \
  (LC08_229082_20200114_toar_B5_hpf + LC08_229082_20200114_toar_B6_hpf) * 1.0"
# Set the color palette
r.colors map=LC08_229082_20200114_NDWI color=ndwi

Mostrar los mapas

ndi = gj.InteractiveMap()
ndi.add_raster()
ndi.show()

NDVI y NDWI a partir de datos Landsat 8

Tarea

Estimar NDVI y NDWI para la otra escena usando el módulo i.vi

Clasificación No Supervisada

• Agrupar las bandas (i.e., hacer un stack): i.group
• Generar firmas para n número de clases: i.cluster
• Clasificar: i.maxlik

Listar los mapas usando un patrón

# list the bands needed for classification
g.list type=raster mapset=. pattern=*_toar*_hpf

Crear un grupo de imágenes o stack

# add maps to an imagery group for easier management
i.group group=l8 subgroup=l8 \
 input=`g.list type=raster mapset=. pattern=*_toar*_hpf sep=","`

Obtener estadísticos -firmas- para las n clases de interés con una muestra de pixeles

# statistics for unsupervised classification
i.cluster group=l8 subgroup=l8 \
 sig=l8_hpf \
 classes=7 \
 separation=0.6

Realizar la clasificación no supervisada de toda la imagen

# Maximum Likelihood unsupervised classification
i.maxlik group=l8 subgroup=l8 \
 sig=l8_hpf \
 output=l8_hpf_class \
 rej=l8_hpf_rej

Mostrar el mapa clasificado

Información derivada adicional

Información derivada adicional podría obtenerse con los siguientes módulos, entre otros:

• medidas de textura: r.texture,
• medidas de diversidad: r.diversity,
• estadísticas locales con información de contexto: r.neighbors,
• transformación tasseled cap: i.tasscap,
• etc.

Clasificación en GRASS GIS

• Topic classification en los manuales de GRASS GIS
• Image classification en la wiki
• Ejemplos de clasificación en el curso dictado en Noruega en 2018
• Clasificación con Random Forest en la presentación del OpenGeoHub Summer School 2018 en Praga
• Detección de cambios con Landsat

Semantic labels

Un concepto bastante nuevo en GRASS GIS son las etiquetas semánticas o semantic labels. Éstas son especialmente relevantes para las imágenes de satélite, ya que nos permiten identificar a qué sensor y banda corresponde una trama determinada. Estas etiquetas son especialmente relevantes a la hora de trabajar con colecciones de imágenes de satélite y también a la hora de clasificar diferentes escenas. Lo veremos más adelante, pero al generar una firma espectral para un determinado conjunto de bandas, puede reutilizarse para clasificar otra escena siempre que las etiquetas semánticas sean las mismas. Cuidado: aunque es posible reutilizar las firmas espectrales para cualquier escena con las mismas bandas, los cambios temporales (estaciones, impacto meteorológico) limitan su aplicabilidad sólo a escenas obtenidas más o menos al mismo tiempo.